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The stellar luminosity and depth of the convective envelope vary rapidly with mass for G- and K-type main sequence stars. 
In order to understand how these properties influence the convective turbulence, differential rotation, and meridional cir- 
culation, we have carried out 3D dynamical simulations of the interiors of rotating main sequence stars, using the anelastic 
spherical harmonic (ASH) code. The stars in our simulations have masses of 0.5, 0.7, 0.9, and 1.1 M©, corresponding to 
spectral types K7 through GO, and rotate at the same angular speed as the sun. We identify several trends of convection 
zone properties with stellar mass, exhibited by the simulations. The convective velocities, temperature contrast between 
up- and downflows, and meridional circulation velocities all increase with stellar luminosity. As a consequence of the 
trend in convective velocity, the Rossby number (at a fixed rotation rate) increases and the convective turnover timescales 
decrease significantly with increasing stellar mass. The 3 lowest mass cases exhibit solar-like differential rotation, in a 
sense that they show a maximum rotation at the equator and minimum at higher latitudes, but the 1.1 M© case exhibits 
anti-solar rotation. At low mass, the meridional circulation is multi-cellular and aligned with the rotation axis; as the mass 
increases, the circulation pattern tends toward a unicellular structure covering each hemisphere in the convection zone. 
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1 Introduction 

The sun and sun-like stars with convection zones in their 
outer envelopes, have long been known to exhibit emission 
line and X-ray activity, associated with hot gas in chro- 



mospheres, transition regions, and coronae (e.g., Pizzolato 



eTaLl [2003| |Strassmeier eTal) [T9901 [Wright et al.| |201 \\ . 
Understanding this activity, and understanding how the so- 
lar case relates to the activity observed accross the HR dia- 
gram is a long-standing puzzle. It is clear that the existence 
of hot gas above the photosphere is related to magnetic pro- 
cesses associated with the convection zone itself. The mag- 
netic field of low-mass main sequence stars is generally be- 
lieved to be generated by dynamo processes, which derive 
their properties from convective motions, differential rota- 
tion, and meridional circulations in the convection zones, as 
well as from the interaction between the convection zone 
and the radiative interior. 

Recent observations, by either spectropolarimetry (e.g., 



Donat i et al.| [2003 ), doppler imaging (e.g., |Barnes et al. 
2005 1), or m onitoring of various activity indicators (e.g., Bal 



iunas et al.l[T995l|Donahue et al.l[T996l|Lovis et all[2DlT] 
Plan et al.| [20091 |Saar & Brandenburg! [T999] ), show that 
solar-like stars possess activity cycles and differential rota- 
tion, analogous to the sun. Solar analogues are even starting 
to be discovered ( [Petit et al.[|2008j ). 

* Corresponding author: e-mail: sean.matt@cea.fr 



In order to gain a better theoretical understanding of 
how convective properties depend upon stellar parameters, 
we carry out 3D numerical dynamical simulations of con- 
vective envelopes of solar-like stars. As a first step, we model 
here the convective regions only, neglecting the effects of 
the interface region with a radiative envelope, and restrict 
ourselves to relatively slow (solar) rotation rates. Specifi- 
cally, we simulate the convection dynamics for 4 main se- 
quence stars with masses of 0.5, 0.7, 0.9, and 1.1 solar masses, 
spanning spectral types GO to K7. This mass interval ex- 
hibits a large range of the physical properties of convective 
envelopes (such as the depth, physical size, mass, and den- 
sity), as well as in the overall stellar luminosity transported 
by convection. These types of stars are also targets for aster- 
oseismic studies (e.g., Verner et al. 201 1| ), which have the 
potential to give precise measurements of stellar properties 
for large numbers of stars. The goal here is to determine how 
the convection, differential rotation, and meridional circula- 
tion is influenced by stellar mass, and to see if general trends 
or scaling laws can be extracted that will guide a deeper un- 
derstanding of the inner hydrodynamics of these stars. The 
present study lays the groundwork for later studies to con- 
sider (e.g.) faster rotation rates, convection-radiation zone 
interface dynamics, and the dynamo generation of magnetic 
fields in stars in this mass range. 

Section[2]contains a description of our simulation method 
and presents a comparison of the overall structures of each 
star. Section [3] describes the main results of our 3D simula- 
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tions, focusing on both the convective turbulence properties, 
as well as the differential rotation and meridional circula- 
tion flows. A summary and brief discussion is contained in 
section [4] 

2 Simulation method 



We use the anelastic spherical harmonic (ASH) code ( Clune 



et al.| |1999| ) to compute the 3 -dimensional and turbulent 
flows in convectively unstable stellar envelopes. This code 
has been extensively tested and used for computing several 



aspects of the solar interior (e.g., Browning et al.| |2006 
Brun et al.l [20041 |Brun & Toomrel |2002l |DeRosa et al. 



20021|Miesch et al.l|2006|), rapidly rotating young stars (|Bal 



lot etal.|[2007l |Brown, 2009; Brown etlH [2008) [2011]), the 



convective cores of massive stars ( Browning et al. 2 004a|b 
Feathersto ne et al.| |2009| ), fully convective low mass stars 
(Browning 2008| ), red giant stars ( [Brun & Palacios 2009| ), 
and pre-main-sequence stars (Bessolaz & Brun, in this vol- 
ume; |Bessolaz & Brun||2011| ). We briefly describe the ba- 
sic aspects of the code here, but the reader can find further 
details of the code in those previous works (see especially, 
|Brun et al.||2004||Clune et al.||l999) . 

The code solves the fluid equations, under the anelas- 
tic approximation, in a computational domain consisting of 
a spherical shell and in a rotating reference frame. Under 
the anelastic approximation, sound waves are filtered out 
and assumed to have a negligible effect on the dynamics, 
in order to allow for a larger computational timestep. This 
approximation is appropriate (e.g.) in the interiors of stars 
because typical motions are highly sub-sonic. The thermo- 
dynamic variables are linearized with respect to a spheri- 
cally symmetric background state with a density p, pressure 
P, temperature T, and specific entropy S and fluctuations 
about the background state of p, P, T, and S. The time- 
dependent equations describe the conservation of mass, mo- 
mentum, and entropy expressed as 



V • (pv) = 0, 



(1) 



' 'dv \ 
p[ — + (v-V)v + 2Q xvj = 

-VP + pg - [VP - pg] - V • X>, (2) 



9 T 



as 
at 



■pTv-V(S + S) 



+ V • \n r pc v V(T + T) + pf^VS + k VS)] 
+ 2pu [eijdj - 1/3(V • vf] , (3) 

where v = (v r , vq, v^) is the velocity in the rotating frame 
in spherical coordinates, tt = ft e z is the angular rota- 
tion rate of the reference frame, g is the acceleration due to 
gravity, K r is the radiative diffusivity, and c p is the specific 
heat at constant pressure. The term T> is the viscous stress 
tensor, with the components 

V l3 = -2pv[e l3 - 1/3(V • v)5 ij ], (4) 



where e i3 - is the strain rate tensor, and v, k, and are effec 
tive eddy diffusivities. The code also uses a linerized equa 
tion of state, 

p _ P T _ P S 
p~P~T~^P~7 p ' 
and the ideal gas law, 



(5) 



P = KpT, 



(6) 



where 7 is the ratio of specific heats (we use 7 = 5/3), and 
1Z is the gas constant. For all 4 stars in our study, the lumi- 
nosity has reached a value of more than 99.9% of the total 
stellar luminosity by the base of the convection zone, so we 
do not need to include any energy generation by nuclear 
burning within our computational domain. 

We set up four different models in which the domain 
boundaries and stratification coincides with ID models of 
the convection zones of main sequence stars with differ- 
describes the global properties of 



ent masses. Section 2.1 



these stars, and sections 



272] and [23] describe the initial and 
boundary conditions used in our ASH models, as well as the 
method for evolving the simulations to a fully-convective, 
statistical steady- state. 

2.1 ID stellar structure 

In order to define the background structure for our 3D mod- 
els, we use the ID stellar evolution code CESAM ( |Morel[ 
|1997| ). With CESAM, we computed the evolution of four 
stars with masses of 0.5, 0.7, 0.9, and 1.1 M until the age 
of 4.6 Gyr. This age is approximately equal to that of the 
sun, so our cases can be compared with the many previous 
results of solar studies. Also, this age is appropriate for the 
slow (solar) rotation rates considered here. For all four stars, 
we assumed the same initial metallicity of (X, Y, Z) « (0.71, 
0.27, 0.02) and a mixing length parameter of 1.77, chosen 
to best represent the solar case. Table [T] lists the key global 
properties of our four stars, at the age of 4.6 Gyr, which are 
also graphically represented in Figures [T}|3] 

Figure [T] shows the total luminosity (solid line), as well 
as the effective temperature (dashed line) as a function of 
mass, for stars in the mass range of our models. The "X" 
symbols mark the values for the 4 stars modeled here, which 
are also listed in table [T] Also indicated on the plot are the 
approximate spectral types, corresponding to the effective 
temperature. The Figure demonstrates the steep rise in lu- 
monosity with mass expected for main sequence stars (the 
stars here approximately follow P* ex M^ 6 ), such that the 
1.1 Mq star is 40 times more luminous than the 0.5 M 
star. 

Figure [2] shows the photospheric radius P* (upper solid 
line) and the radial location of the base of the convection 
zone R cz (lower solid line) as a function of mass. The stel- 
lar radius is strong function of stellar mass in the plotted 
range, with the 1.1 M star being 2.8 times bigger than 
the 0.5 Mq star. The radial extent of the convection zones 
for these stars is represented by the region between the two 
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Table 1 Global properties of the 4 stars used in our ASH models, computed with the CESAM stellar evolution code and 
at an age of 4.6 Gyr. We adopt M© = 1.989 x 10 33 g, R Q = 6.9599 x 10 10 cm, and L© = 3.846 x 10 33 erg s -1 . 
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Fig. 1 Stellar luminosity in solar units (solid line, left 
scale) and effective temperature (dashed line, right scale) 
as a function of stellar mass, for the mass range spanned 
by our models, and as computed using the CESAM stellar 
evolution code. The "X" symbols indicate the values for the 
4 stars used in our 3D ASH simulations. The approximate 
spectral types are also indicated next to each of these. 



0.8 

M*/M sun 

Fig. 2 Stellar radius (upper line) and radial location of the 
base of the convection zone (lower line) in solar units and 
as a function of stellar mass, as computed by the CESAM 
stellar evolution code. The convection zone exists between 
these two lines, while a convectively stable radiation zone 
lies below. The "X" symbols indicate the upper and lower 
radial boundaries of the computational domain used in our 
ASH simulations, for the 4 cases considered. 



solid lines, and the region below contains a convectively sta- 
ble radiation zone. It is clear that, in the mass range shown, 
the thickness of the convection zone increases slightly with 
increasing stellar mass. However, since the stellar radius in- 
creases more rapidly, the fractional size of the convection 
zone decreases with increasing mass. Thus, the convection 
zone thickness ranges from 44% to 25% of the stellar radius, 
for the 0.5 and 1.1 M stars, respectively (for the sun, this 
value is approximately 30%). The "X" symbols indicate the 
upper and lower boundaries of our 3D simulation domains 
for the 4 stars simulated (discussed below). 

Figure[3]shows the extent of convection zones, expressed 
in mass coordinates and normalized to the mass of each star, 
as a function of stellar mass. The lines show the total mass 
enclosed by the stellar surface (upper line) and enclosed by 
the location of the base of the convective envelope (lower 
line). The amount of mass contained in the convection zone, 
M cz , is the difference between these two lines and is listed 
in Table Q] It is clear that the convection zone mass is a 



strong function of stellar mass in the range considered. Stars 
with slightly lower mass will be fully convective, while sig- 
nificantly more massive stars will not have convective en- 
velopes. The convection zone mass varies from 36% to 1% 
of the stellar mass, for the 0.5 and 1.1 M stars, respectively 
(the solar value is approximately 3%). The "X" symbols in- 
dicate the simulation domain boundaries, expressed in mass 
coordinates, for the 4 stars simulated (discussed below). 

The last two columns of table [T] show the temperature 
and the mass density at the base of the convection zone for 
each case. Although it is true that the temperature and den- 
sity increases with mass at the stellar center (R = 0), these 
quantities decrease with increasing mass at the location of 
R cz . In other words, for increasing stellar mass, the con- 
vection zone generally occupies a more tenuous and cooler 
outer layer of the star. 
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Table 2 Simulation parameters for each case. The nuber of radial, latitudinal, and longitudinal gridpoints are N r , N$, 
and N^. The outermost radius of the simulated spherical shell is R ou t, and H p (R out ) is the density scale height there. The 
radial size of the domain is L = R out — R cz . The viscosity at R out is z^top and it varies in the domain with the inverse 
square root of the background density. All stars have a Prandtl number P r = v/k = 0.25 and rotate at the solar rate, 
= 2.6 x 10 -6 rad/s, corresponding to a rotation period of 2tt/Qo & 28 days. Also listed are the Rayleigh number R a = 
(-dp/dS)ASgL 3 /(pvK), the Taylor number T a = 4^ 2 L 4 /^ 2 , and the convective Rossby number R oc = (R a /T a P r ) 1/2 , 
all evaluated at the midlevel of the domain. 
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Fig. 3 Mass contained beneath the stellar surface (upper 
line) and location of the base of the convection zone (lower 
line) as a function of stellar mass, as computed by the CE- 
SAM stellar evolution code. The convection zone exists be- 
tween these two lines, while a convectively stable radiation 
zone lies below. The "X" symbols indicate the upper and 
lower radial boundaries of the computational domain used 
in our ASH simulations, for the 4 cases considered. 



2.2 Initial and boundary conditions for 3D 
simulations 

We use the ID stellar structure models described above for 
the spherically symmetric, initial conditions of the the 3D 
ASH models. Table[2]lists the key properties of each of the 4 
simulations. In the ASH models presented here, the domain 
consists of the stellar convection zone only. Thus, the lower 
radial boundary of our simulated spherical shells coincides 
with the base of the convection zone, R cz , for each star (see 
table [T] and the lower row of "X" symbols in Figs. [2] and 
[3}. We chose the location of the outer boundary such that 
the mass density there is a factor of 100 smaller than at the 
base of the convection zone. This places the outer boundary 



somewhat below the photosphere of the star. The location 
of the outer boundary ,R ou t is listed in table [2|and shown as 
the upper row of "X" symbols in and Figures[2]and[3] 

In each of our 4 cases, the gravity (g) is taken directly 
from the corresponding CESAM model, and the entropy 
gradient (dS/dr) is initialized with a constant value equal 
the (superadiabatic) entropy gradient in the middle of the 
convection zone in the CESAM model. The background 
thermodynamic variables (p, P, T, and S) are set accord- 
ing to this entropy gradient and to be in hydrostatic bal- 
ance with the stellar gravity, while the fluctuating thermody- 
namic variables (p, P, T, and S) are initially zero. The ini- 
tial entropy gradient and choice of diffusivities (see below) 
ensures that the initial state is convectively unstable and has 
a Rayleigh number near or above the critical Rayleigh num- 
ber necessary for convection. The velocity in the rotating 
frame (v) is given an initial random perturbation, with negli- 
gible kinetic energy, in order to initiate convective motions. 

For both the inner and outer boundaries, we use stress- 
free and impenetrable boundary conditions on the velocity 
and hold the entropy gradient fixed at the initial value. The 
radiative diffusivity Kjf IS chosen so that the ASH model has 
the same radiative flux at all radii as the corresponding CE- 
SAM model. The initial stellar structure and boundary con- 
ditions ensure that the energy flux into the domain at the 
bottom of the convection zone (given entirely by the radia- 
tive flux) is constant in time. At the top of the domain, the 
density scale height becomes small (see table [2]), and since 
this approximately determines the convection cell size scale, 
it becomes numerically challenging to resolve the convec- 
tive motions there. Furthermore, the impenetrable bound- 
ary condition precludes any convective enthalpy flux from 
escaping the top of the domain. To address both of these 
issues we introduce a diffusive energy flux (the term pro- 
portional to fto in eq. (3)) that is assumed to represent an 
unresolved and spherically symmetric enthalpy flux carried 
by small-scale convection near the top of the domain (here- 
after "unresolved eddy flux"). In all 4 cases, we chose so 
that this flux is negligible in the bulk of the convection zone, 
but increases as a smooth function near the outer boundary, 
such that the flux leaving the domain equals the stellar en- 
ergy flux and is constant in time. 
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Fig. 4 Energy flow in 0.5 (left) and 1.1 M (right) cases as a function of radius in the computational domain. Shown are 
the spherically- and time-averaged luminosities, after the simulations have reached a statistical steady-state. 



In each case, we chose the viscosity at the outer part of 
the domain, is top , in order to ensure that the resulting con- 
vective turbulence has a significant Reynolds number. As 
shown in section[3] a higher stellar luminosity leads to larger 
convective velocities, and thus a larger viscosity is needed 
to keep the Reynolds numbers comparable (explaining the 
trend of is t0 p apparent in table [2]). Within the domain, the 
viscosity in each case depends only upon radius, and it var- 
ries as the inverse square root of the background density. 
Scaling the viscosity in this way ensures a significant tur- 
bulence level at all radii and is used in most previous ASH 
models in the literature. The thermal diffusivity n equals 4u 
everywhere in the domain, giving a constant Prandtl num- 
ber of 0.25 for all cases. The rotation rate of the reference 
frame (£lo), which equals the rotation rate of the star, is set 
to the solar rate for all 4 stars. 

The Rayleigh numbers (calculated using the steady- state 
value of the entropy in the simulations), Taylor numbers, 
and convective Rossby numbers for each case are also given 
in table [2] The Rayleigh and Taylor numbers, show a strong 
decrease with increasing stellar mass, which is primarily 
explained by the trend in the diffusivities. If the diffusivi- 
ties were constant, these numbers would increase with mass 
for these stars. The convective Rossby number generally in- 
creases with mass, meaning that for the same rotation rate, 
the coriolis force has less influence on the dynamics. 

2.3 Reaching a statistical steady-state and energy flux 
balance 

At the start of the simulation, the star is in a quiescent state. 
Thus, there is no significant enthalpy flux, and the system is 
not initially in radial energy flux balance. In this situation, 
the evolution of the system, according to equations ([T])-([3]), 
ensures that the entropy gradient evolves toward an energy 



flux balanced state. Furthermore, since the gas in the com- 
putational domain is unstably stratified (dS/dr < 0), with 
a large Rayleigh number (see table [2]), significant convec- 
tive motions begin rapidly after the start of the simulations. 
Once vigorous convection begins, an energy flux balance is 
achieved within a few convective turnover times (months, 
typically) and maintained for the duration of the simula- 
tions. In this state, the net energy transport across the do- 
main is constant at all radii, when averaged over several 
convective turnover times. 

Figure [4] shows the spherically- and time-averaged lu- 
minosity, as a function of radius in the whole domain, for 
the highest and lowest mass stars in this study. These are 
shown after the simulated stars have evolved for several 
years from the initial state. The boundary conditions ensure 
a fixed radiative energy flux into the domain at the inner 
boundary and a fixed unresolved eddy flux out of the do- 
main at the outer boundary, both with luminosities equal 
to that of the modeled stellar luminosity. Within the do- 
main, the radiative energy flux is nonzero, but it decreases 
in importance with increasing radius in the convection zone. 
The enthalpy flux, associated with convective motions, car- 
ries the bulk of the remaining energy flux from the lower 
boundary, across the convection zone, to the outer bound- 
ary. Due to the asymmetry between broad, slow upflows 
and narrow, fast downflows, there is generally a net neg- 
ative (downward) flux of kinetic energy in the convection 
zone (e.g., as in the right panel of Fig. |4j Catta neo et al. 
|1991||Hurlburt et al.| [T986 ). As a result, the steady-state en- 
thalpy luminosity (dash-triple-dotted line) peaks at a value 
significantly larger than L*, in order that the net energy flow 
across the convection zone is constant and equal to the stel- 
lar luminosity. As demonstrated by the two extreme cases 
in Figure [4] the presence of a super-luminal enthalpy flux is 
increasingly important with mass, in our 4 simulations. It is 
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Fig. 5 Integrated kinetic energies in the domain, divided 
by the total volume, as a function of time, for the 0.5 M 
case. Shown are the kinetic energies associated with the dif- 
ferential rotation (DRKE), meridional circulation (MCKE), 
convective motions (CKE), and the total (solid line). 



not clear whether this is due to an intrinsic property of the 
stars (such as the luminosity), or whether it only depends 
upon (e.g.) the viscosity, which is systematically different 
in each of our cases (table [2}. However, it does seem that 
the faster flows found in the more massive stars (discussed 
below) yields a larger negative kinetic energy flux, since the 
difference between the broad- slow upflows and narrow-fast 
downflows is accentuated. 

Although the energy flux balance is established rela- 
tively rapidly in the convection zone, the statistical prop- 
erties of the turbulence can take longer to settle down, as 
they evolve in response to (e.g.) Reynolds stresses and the 
dissipation of energy on small scales. In addition, these sys- 
tems typically exhibit global, axisymmetric flows such as 
differential rotation and meridional circulation. These flows 
are also influenced by the Reynolds stresses, and the merid- 
ional circulation typically has velocities smaller than the 
convective velocities. Thus, the global flows settle down on 
timescales of the order of 10 times the convective turnover 
time. In order to ensure the simulations have time to estab- 
lish these global flows and are near their turbulent statistical 
steady-state, we have run each model for a timescale com- 
parable to the global viscous dissipation timescale (L 2 /v). 

As an example, Figure [5] shows the globally integrated 
kinetic energies in the 0.5 M case, as a function of time 
from the start of the simulation. The evolution during the 
first year is characterized by an initially linear (exponen- 
tial) growth of the convective instability, followed by a non- 
linear saturation phase, then by relaxation oscillations until 
it settles down into a statistically steady state. The longer 
term evolution exhibits a growth in the differential rotation, 



which also reaches a steady state. It is clear from Figure [5] 
that there is little change during the last 5000 days of evolu- 
tion. The results in Figure [4] and presented below are shown 
during this mature state of the systems. 



3 3D simulations 

The ending time of each of the 4 simulations (£/) is listed in 
table [3] In this section, we present the convective properties 
and the global flows existing in each case. 



3.1 Properties of convection 

Figure [6] shows the radial velocity on slices of constant ra- 
dius at slightly below the top of the domain and at the mid- 
dle of the domain, for the 0.5 (top), 0.7 (middle), and 1.1 
(bottom) Mq cases. Also shown are the slices of tempera- 
ture at the middle of the domain (right panels). The temper- 
ature slices generally appear less structured than the radial 
velocity slices because the Prandtl number is less than unity 
(0.25), so that the thermal difussivity (k) is larger than the 
kinematic viscosity (y). 

The convective patterns and characteristic eddy sizes are 
influenced strongly by the density scale height, as well as by 
the contrast between horizontal and vertical velocity (see, 
e.g., Bessolaz & Brun 201 1) . Near the surface, where the 
density scale height is smallest, the shell slices exhibit typ- 
ical convective patterns characterized by hot, broad/patchy 
upflows surrounded by cool, narrow downflow lanes (e.g., 
|Cattaneo etaE) [1991] |Miesch et aT) [2008] ). Deeper in the 
convection zone, the downflow lanes often merge and gen- 
erally lose their connectivity with respect to the surface pat- 
terns, and the convection has a less patchy appearance. The 
near- surface shell slices (left panels of Fig. [6]) are shown at a 
location where the unresolved eddy flux dominates the other 
fluxes (see Fig.|4j. The unresolved eddy flux is a spherically 
symmetric quantity, and it does not affect the shape of the 
convective patterns, but it does influence the level of (re- 
solved) turbulence, since the unresolved eddy flux reduces 
the convective driving. However, at the location of the near- 
surface shell slices shown, all cases are still turbulent, with 
a significant convective (enthalpy) flux and Reynolds num- 
bers in the range of 10 — 30. 

In all of our cases, the convective patterns evolve on typ- 
ical timescales of a few weeks to a few months, with cells 
merging, splitting, and disappearing. The evolution of the 
convective motions in the 1 . 1 M case generally evolve on 
timescales of a factor of a few times shorter than in the 0.5 
Mq case. The convective patterns are partly advected by 
the local shear and, depending on whether the flow is pro- 
grade of retrograde, are tilted to the right or the left near the 
equator. At higher latitude, beyond the tangent cylinder (an 
imaginary circle crossing the upper surface with a cylindri- 
cal radius of R cz ), the convective cells are less aligned with 
the rotation axis and exhibit a more patchy behavior. 
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Table 3 Simulation results for each case. The amount of evolution time since the start of the simulation is denoted tf. 
Temporal averages of the rms components of velocity v r ,v^,VQ, speed v, and fluctuating velocities and v' are evaluated 
at the midlevel of the domain and given in units of m s _1 . The temporally averaged rms temperature fluctuation at midlevel 
is T. Also listed are the rms Reynolds number R e = v'L/v, Rossby number R Q = v'/(2QoL), and Peclet number 
P e = v'L/k, evaluated at the midlevel, and the convective turnover timescale r to = L/v r . The is the difference in 
angular rotation rate at the outer boundary, between latitudes of 0° and 60°. The AT is the difference in the temporally 
and azimuthally averaged temperature at the base of the convection zone (R cz ), between latitudes of 60° and 0°. Finally, 
v mc is the rms meridional circulation speed, calculated from temporally and azimuthally averaged poloidal velocity and by 
taking the rms of all values at a constant midlevel radius. 
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Fig. 6 Slices at constant radius (Mollweide view) for the 0.5 M case (top), 0.7 M case (middle), and 1.1 M case 
(bottom), showing the radial velocity near the surface (left panels) and in the midlevel (middle panels) of the domain and the 
temperature at midlevel (right panels). In the temperature slices, the azimuthal average has been subtracted, to emphasize 
the temperature fluctuations. The dotted line corresponds to the stellar surface and the dashed line to the equator. 
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Fig. 7 Azimuthally and temporally averaged Rotation frequency in a meridional slice through the domain (contour/color 
plots; rotation axis is vertical) and shown along radial lines at different latitudes (line plots, latitude given in degrees). 
Shown are the 0.5 (top left), 0.7 (top right), 0.9 (bottom left), and 1.1 (bottom right) M cases. 



As evident in the Figure, the 0.5 M case exhibits con- 
vective patterns near the equatorial region that are inhomo- 
geneous, consisting of vigorous convection over some range 
of longitudes and relatively quiescent regions at other lon- 
gitudes. Such "active nest" convective patterns have been 
discussed by |B allot et all] p007] ) and |Brown et aL] ( [2008] ) 
and are likely related to our choice of low Prandtl number. 
The remaining 3 cases, exhibit vigorous convection within 
the entire volume of the convection zone, as seen for the 0.7 
and 1 . 1 M cases shown in the bottom two rows of Figure 

m 

Table |3] lists the rms velocities at the middle of the com- 
putational domain for each case. There is a clear trend of in- 
creasing convective velocities with increasing stellar mass. 
As evident in the right panels of Figure [6] and in the values 
of the rms temperature (T) listed in table [5J there is also a 
strong trend of increasing temperature contrast between up- 



and downflows, with increasing mass. This trend is steeper 
if one considers the temperature fluctuation relative to the 
mean temperature, since the mean background temperature 
in the 1.1 M convection zone midlevel is ~3 times colder 
than for the 0.5 M star (see jp.l|). Furthermore, the aver- 



age mass density in the convective envelope decreases with 
increasing stellar mass (the midlevel density in the 1.1 M 
case is 280 times lower than in the 0.5 M case), which 
decreases the efficiency of convective energy transport. In 
spite of the lower density, the convection in higher mass 
stars is able to carry a larger luminosity than the lower mass 
stars, due to a combination of larger stellar radii, higher con- 
vective velocities, and higher temperature contrasts between 
up- and downflows. 

Table[3]also lists the Reynolds, Rossby, and Peclet num- 
bers for each model. The Reynolds (and Peclet) have a range 
of a factor of a few, with no trend in mass. These values pri- 
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marily reflect our choice of u top for each case (tab. [2]) rather 
than any intrinsic property of the stars. For our choices of 
z^top, the 0.5 Mq case has the highest Reynolds number. On 
the other hand, the Rossby numbers do not depend strongly 
on our choice of diffusivities, and they show a strong de- 
pendence on mass, primarily due to the trend in convective 
velocity with mass, at a constant rotation period. Finally, ta- 
ble lists a characteristic convective turnover timescale for 
each case. This also shows a strong trend with mass, primar- 
ily due to the trend in convective velocities. 

3.2 Differential rotation and meridional circulation 

Figure [7] shows the differential rotation in all 4 cases. The 
2D plots in the Figure exhibit isorotation contours that are 
nearly aligned on cylinders in all cases. This behavior is 
typical for simulations such as these that have no enhanced 
latitudinal entropy gradient present at the lower boundary 
of the convection zone ( [Ballot et al.| |2007| |Miesch et al.| 
[2006] ). 

The angular rotation rate in the 0.5 M case has a min- 
imum value at mid latitudes, a local maximum in the outer 
equatorial region, and a global maximum near the pole. A 
similarly "banded" differential rotation pattern is evident in 
a number of previously published ASH simulations (e.g. 



Bessolaz & Brun 20TT] Browning 2008). The differential 



rotation in the 0.7 and 0.9 M cases are the most solar-like, 
with the fastest angular rotation rate at the equator and slow- 
est at higher latitudes. However, the 1.1 M case exhibits 
anti-solar rotation, where the slowest angular velocity is at 
the equator. A reversal of the sense of differential rotation 
was also observed in the simulations of Bessolaz & Brunl 
( |2011| ) to be an effect of the thickness of the convection 
zone. In that study, the stellar structure and luminosity was 
held fixed, while only the thickness of the convection zone 
was varried. The thickness of the convection zone in our 1.1 
Mq case (as a fraction of R*) is comparable to the models 



of Bessolaz & Brun ( 201 1 ) showing anti-solar rotation. This 
would suggest that the switch from solar to anti-solar may 
have more to do with convection zone thickness than with 
stellar structural properties. At the same time, our choice of 
parameters may also have influenced this outcome. In par- 
ticular, the 1.1 M case has the largest Rossby numbers 
(see tabs. [2] and [3}, which means that the convection is less 
influenced by rotation than in the other cases. 

Each case exhibits a time-averaged, axisymmetric, lati- 
tudinal temperature gradient. For the 3 highest mass cases, 
this is characterized by a monotonic change from equator 
to pole. The 1.1 M case has a hotter equator, while the 
0.7 and 0.9 M cases have hotter poles (as in the sun). The 
0.5 Mq case has a more complex temperature pattern with a 
maximum at mid latitudes, while the latitudinal entropy gra- 
dient in all cases follows a more monotonic behavior from 
equator to pole. The second-to-last column of table [3] shows 
the temporally and azimuthally averaged temperature differ- 
ence between the equator and a latitude of 60° , at the base 
of the convection zone. It is evident that the temperature 




max 



mm 



Fig. 8 Contours of the stream function of the meridional 
flow for the 0.5, 0.7, 0.9, and 1.1 M cases (left to right), 
where the physical sizes of the stars are shown to scale. The 
contours follow streamlines in the meridional plane from 
temporally and azimuthally averaged data. The colors indi- 
cate clockwise (blue) and counter-clockwise (red) circula- 
tion, and the stellar rotation axis is vertical. 



gradients are steaper for higher mass stars. The more com- 
plex temperature structure of the lowest mass case is con- 
sistent with the banded differential rotation pattern, and the 
reversed polarity of the temperature variation in the highest 
mass case is consistent with anti-solar differential rotation. 

The fractional differential rotation of each case, mea- 
sured at the top of the domain between latitudes and 60 
degrees, is listed in the final column of table [3] These val- 
ues of Aft/fto in the table, show that the magnitude of the 
differential rotation ranges from 30 to 42% for all cases. 

Figure[8]shows the meridional circulation for all 4 cases. 
The plots are shown to scale, according to the physical size 
of each of the stars. For the 2 highest mass stars, the cir- 
culation pattern is dominated by a global, unicellular flow 
pattern from equator to pole near the surface and from pole 
to equator near the base of the convection zone. However, 
as evident in the Figure, there is a trend for this flow to 
break up into a multi-cellular network of circulation pat- 
terns roughly aligned on cylinders, for the lower mass stars. 
The tendancy for the meridional flow to breakup in this way 
for increasing convection zone depth was also seen in the 
simulations of |Bessolaz & Brun] ( |201 1| ), suggesting that this 
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effect is mostly due to the size of the convection zone, rather 
than with other stellar properties. 

The last column of table|3]shows the characteristic merid- 
ional circulation flow speed evaluated at the midlevel of the 
domain for each case. As with the convective velocities, 
there is a clear trend of increasing meridional circulation 
flow speed with increasing stellar mass. A similar trend is 



also seen in 2D models (e.g., Kuker & Rudiger , 2008 ). Since 
the physical size of the circulating patterns (evident in Fig. 
[8} also increases strongly with mass, this suggests that the 
characteristic timescale for these circulations is not strongly 
dependent on stellar mass. The values listed in table [3] also 
suggest that the ratio of the meridional circulation speed to 
the convection speed (e.g., v mc /v) range from 2% for the 3 
lowest mass stars to 9% for the 1.1 M case. 

4 Summary & discussion 

We have presented results from 3D dynamical simulations 
of the convection zones of 4 sun-like stars, spanning a mass 
range from 0.5 to 1.1 M and rotating at the solar rate. As 



presented in section |2.1| these stars cover a range in lumi- 
nosities of a factor of 40 and have convection zone masses 
ranging from 18% to 1.1% of a solar mass. The physical 
sizes of the convection zones increase slightly with mass 
(ranging from 19% to 31% of a solar radius), although the 
fractional size decreases with mass (ranging from 44% to 
24% of a stellar radius). 

We found that the convective velocities and temperature 
contrast between up- and downflows is a strong function 
of mass. This is due to the strong luminosity dependence 
on mass, coupled with the fact that higher mass stars have 
less dense convective envelopes and thus require more vig- 
orous convection. As a consequence of the trend in convec- 
tive velocities, the Rossby number (for a constant Qo) in- 
creases and the convective turnover time decreases signifi- 
cantly with increasing mass. 

The convective pattern in the 0.5 M case exhibits "ac- 
tive nests" on a finite range of longitudes in the equatorial 
region. The presence of these active nests appears to be re- 
lated to both the Rossby and Prandtl numbers, in a sense that 
they have been seen in simulations with low R Q number and 
with P r = 0.25 ( [Ballot et al.||2007l|Brown et al.||2008] ) but 
not in cases with P r = 1.0 ( [Ballot et alT] [20071 |Bessolaz| 
& Brun, 201 1 1 . This may be due to the relative amplitude 
of the local shear, which tends to be stronger for lower P r . 
Stronger shear disrupts convective motions and can locally 
inhibit radial motions. 

The angular rotation rate has a peak value at the sur- 
face and equator for the 0.5, 0.7, and 0.9 M cases. How- 
ever, our simulations exhibit a differential rotation that is 
"banded," with a minimum at mid latitudes and another max- 
imum at the poles, for the 0.5 M case. The 0.7 and 0.9 M 
cases exhibit the most solar-like differential rotation, with a 
maximum at the equator and minimum near the poles, but 
the 1 . 1 M case exhibits anti-solar differential rotation. The 



magnitude of the differential rotation is comparable in all 
cases. 

The 4 simulations presented here exhibit isorotation con- 
tours that are nearly aligned on cylinders, approximately 
following the Taylor-Proudman constraint. From simulations 
of the solar convection zone (e.g., Miesch et al. 2006 ), we 
know that the Taylor-Proudman constraint can be broken by 
a thermal wind driven by very small latitudinal gradients 
of temperature or entropy, with differences of a few Kelvin 
sufficing to move cylindrical profiles to radially-aligned ro- 
tation profiles. Such gradients and thermal winds do occur 
in the simulations presented here, but the presence of the do- 
main boundary at the base of the convection zone somewhat 
inhibits this effect. By moving the lower domain bound- 
ary deep into the radiation zone,|Brun et al.|p011|) demon- 



strated that the presence of a dynamically self-consistent 
tachocline alleviates the effects of the lower boundary on 
the convection zone dynamics and results in a stronger ther- 
mal wind, isorotation contours more aligned to the radial di- 
rection, and more realistic meridional circulations and con- 
vective penetration. In future work, similar models will be 
developed for the stars in the present study. 

The meridional circulation is dominated by a single, global 
circulation pattern in each hemisphere for the most massive 
cases, but breaks up into smaller scale patterns, aligned with 
the cylindrical z-direction, for lower mass stars. 

Future work with these models should explore cases that 
include, for example, (a) a part of the radiation zone in the 
computational domain, in order to determine how the pres- 
ence of a tachocline affects the convection zone properties, 
(b) the effects of different rotation rates, in particular faster 
rotation, to understand how the properties of these stars dif- 
fer at younger ages, and (c) magnetic fields, in order to ad- 
dress the magnetic dynamo properties of these stars. 
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